####
# Replication code for "Households with solar installations are ideologically diverse and more politically active than their neighbors"
# Mildenberger, Matto; Howe, Peter D.; Miljanich, Chris
# Accepted October 2019, Nature Energy
#### 

library(data.table)
library(ggplot2)
library(stargazer)

# Load "hhdata.Rda" (file contains all solar and control households, including those with missing targetsmart data)
# Load "hhdata.matches.Rda" (file contains only matched pairs of households where both have targetsmart data)

#descriptives for all adopters
stargazer(hhdata[hhdata$type=="solar",c(4,7,8,21,17,18,19,24,20,9,10,11,12,16,15)], summary.stat=c("n", "mean", "sd", "min", "max"), type="latex", header=FALSE)
#descriptives for all matched controls
stargazer(hhdata[hhdata$type=="control",c(4,7,8,21,17,18,19,24,20,9,10,11,12,16,15)], summary.stat=c("n", "mean", "sd", "min", "max"), type="latex", header=FALSE)

#descriptives by quintile
scores.bystrat.type <- hhdata[ , .(.N, mean(tsmart_partisan_score_hhmean, na.rm=TRUE), mean(dem.prop, na.rm=TRUE), mean(rep.prop, na.rm=TRUE), mean(voteg.prop, na.rm=TRUE), mean(votep.prop, na.rm=TRUE)),keyby=.(strat, type)]
names(scores.bystrat.type) <- c("Solar density quintile", "Type", "N", "Partisan score mean", "Prop. Dem.", "Prop. Rep.", "Vote general", "Vote primary")
stargazer(scores.bystrat.type, summary=FALSE, header=FALSE,type="latex")

#number of individuals in each group
sum(hhdata.matches$n.hh[hhdata.matches$type=="solar"])
sum(hhdata.matches$n.hh[hhdata.matches$type=="control"])

## Descriptives for matched file
#descriptives for adopters
stargazer(hhdata.matches[hhdata.matches$type=="solar",c(4,7,8,21,17,18,19,24,20,9,10,11,12,16,15)], summary.stat=c("n", "mean", "sd", "min", "max"), type="latex", header=FALSE)
#descriptives for matched controls
stargazer(hhdata.matches[hhdata.matches$type=="control",c(4,7,8,21,17,18,19,24,20,9,10,11,12,16,15)], summary.stat=c("n", "mean", "sd", "min", "max"), type="latex", header=FALSE)

###
# Created separate paired dataset with only pairs where party registration data are available for both
###
list.s <- hhdata.matches$matchid[is.na(hhdata.matches$dem.prop)==FALSE & hhdata.matches$type=="solar"]
list.c <- hhdata.matches$matchid[is.na(hhdata.matches$dem.prop)==FALSE & hhdata.matches$type=="control"]
hhdata.matches.party <- hhdata.matches[hhdata.matches$matchid %in% list.s,]
hhdata.matches.party <- hhdata.matches.party[hhdata.matches.party$matchid %in% list.c,]

#### 
#tests and figures
####

### Figure 2a in paper
quartz(width=4, height=4) 
p <- ggplot(data=hhdata.matches, aes(x=tsmart_partisan_score_hhmean)) + geom_line(stat="density",aes(color=type)) + theme_bw() + scale_color_brewer(palette="Dark2", name="Household\ntype") + xlab("Partisan score") + scale_x_continuous(breaks=c(0,30,70,100), labels=c("Strong\nRep.", "Lean\nRep.", "Lean\nDem.", "Strong\nDem.")) + theme(panel.grid=element_blank(), legend.position="bottom")
p

ks <- ks.test(hhdata.matches$tsmart_partisan_score_hhmean[hhdata.matches$type=="solar"], hhdata.matches$tsmart_partisan_score_hhmean[hhdata.matches$type=="control"])
tt <- t.test(tsmart_partisan_score_hhmean ~ type, na.rm=TRUE, data=hhdata.matches, paired=TRUE) 



### Figure 2b in paper (party reg)
## combined party registration
temp <- hhdata.matches.party[,.(mean(rep.prop, na.rm=TRUE), sd(rep.prop, na.rm=TRUE)/sqrt(.N), mean(rep.prop, na.rm=TRUE)-2*(sd(rep.prop, na.rm=TRUE)/sqrt(.N)), mean(rep.prop, na.rm=TRUE)+2*(sd(rep.prop, na.rm=TRUE)/sqrt(.N))), by=.(type)]
names(temp) <- c("type", "mean", "se", "lower", "upper")
temp$party <- "Republican"
temp2 <- hhdata.matches.party[,.(mean(dem.prop, na.rm=TRUE), sd(dem.prop, na.rm=TRUE)/sqrt(.N), mean(dem.prop, na.rm=TRUE)-2*(sd(dem.prop, na.rm=TRUE)/sqrt(.N)), mean(dem.prop, na.rm=TRUE)+2*(sd(dem.prop, na.rm=TRUE)/sqrt(.N))), by=.(type)]
names(temp2) <- c("type", "mean", "se", "lower", "upper")
temp2$party <- "Democratic"
temp3 <- hhdata.matches.party[,.(mean(unaff.prop, na.rm=TRUE), sd(unaff.prop, na.rm=TRUE)/sqrt(.N), mean(unaff.prop, na.rm=TRUE)-2*(sd(unaff.prop, na.rm=TRUE)/sqrt(.N)), mean(unaff.prop, na.rm=TRUE)+2*(sd(unaff.prop, na.rm=TRUE)/sqrt(.N))), by=.(type)]
names(temp3) <- c("type", "mean", "se", "lower", "upper")
temp3$party <- "Unaffiliated/\nno party"
temp <- rbind(temp,temp2,temp3)

quartz(width=5, height=3)
p <- ggplot(data=temp, aes(y=mean, x=type)) + geom_point(aes(color=type)) + theme_bw() + geom_errorbar(aes(ymax=upper, ymin=lower, color=type), width=.2, size=.2) + scale_color_brewer(palette="Dark2", name="Household\ntype") + ylab("Proportion \nby party registration") + facet_wrap(~party) + ylim(0,.48) + xlab("Household type") + theme(panel.grid=element_blank())
p



### Figure 5, party registration by solar density stratum
temp <- hhdata.matches.party[,.(mean(rep.prop, na.rm=TRUE), sd(rep.prop, na.rm=TRUE)/sqrt(.N), mean(rep.prop, na.rm=TRUE)-2*(sd(rep.prop, na.rm=TRUE)/sqrt(.N)), mean(rep.prop, na.rm=TRUE)+2*(sd(rep.prop, na.rm=TRUE)/sqrt(.N))), by=.(type, strat)]
names(temp) <- c("type", "strat", "mean", "se", "lower", "upper")
temp$party <- "Republican"
temp2 <- hhdata.matches.party[,.(mean(dem.prop, na.rm=TRUE), sd(dem.prop, na.rm=TRUE)/sqrt(.N), mean(dem.prop, na.rm=TRUE)-2*(sd(dem.prop, na.rm=TRUE)/sqrt(.N)), mean(dem.prop, na.rm=TRUE)+2*(sd(dem.prop, na.rm=TRUE)/sqrt(.N))), by=.(type, strat)]
names(temp2) <- c("type", "strat","mean", "se", "lower", "upper")
temp2$party <- "Democratic"
temp3 <- hhdata.matches.party[,.(mean(unaff.prop, na.rm=TRUE), sd(unaff.prop, na.rm=TRUE)/sqrt(.N), mean(unaff.prop, na.rm=TRUE)-2*(sd(unaff.prop, na.rm=TRUE)/sqrt(.N)), mean(unaff.prop, na.rm=TRUE)+2*(sd(unaff.prop, na.rm=TRUE)/sqrt(.N))), by=.(type, strat)]
names(temp3) <- c("type", "strat", "mean", "se", "lower", "upper")
temp3$party <- "Unaffiliated/\nno party"
temp <- rbind(temp,temp2,temp3)

quartz(width=6, height=3)
p <- ggplot(data=temp, aes(y=mean, x=strat, group=type)) + geom_point(aes(color=type), position=position_dodge(width=0.5)) + theme_bw() + geom_errorbar(aes(ymax=upper, ymin=lower, color=type), width=.2, size=.2, position=position_dodge(width=0.5)) + scale_color_brewer(palette="Dark2", name="Household\ntype") + ylab("Proportion \nby party registration") + facet_wrap(~party) + ylim(0,.75) + xlab("Quintile of census tracts by density of solar installations") + theme(panel.grid=element_blank())
p



### Figure 4 in paper (voting behavior)
#voting history (combined)
temp <- hhdata.matches[,.(mean(voteg.prop, na.rm=TRUE), sd(voteg.prop, na.rm=TRUE)/sqrt(.N), mean(voteg.prop, na.rm=TRUE)-2*(sd(voteg.prop, na.rm=TRUE)/sqrt(.N)), mean(voteg.prop, na.rm=TRUE)+2*(sd(voteg.prop, na.rm=TRUE)/sqrt(.N))), by=.(type)]
names(temp) <- c("type", "mean", "se", "lower", "upper")
temp$election <- "General"
temp2 <- hhdata.matches[,.(mean(votep.prop, na.rm=TRUE), sd(votep.prop, na.rm=TRUE)/sqrt(.N), mean(votep.prop, na.rm=TRUE)-2*(sd(votep.prop, na.rm=TRUE)/sqrt(.N)), mean(votep.prop, na.rm=TRUE)+2*(sd(votep.prop, na.rm=TRUE)/sqrt(.N))), by=.(type)]
names(temp2) <- c("type", "mean", "se", "lower", "upper")
temp2$election <- "Primary"
temp3 <- hhdata.matches[,.(mean(votem.prop, na.rm=TRUE), sd(votem.prop, na.rm=TRUE)/sqrt(.N), mean(votem.prop, na.rm=TRUE)-2*(sd(votem.prop, na.rm=TRUE)/sqrt(.N)), mean(votem.prop, na.rm=TRUE)+2*(sd(votem.prop, na.rm=TRUE)/sqrt(.N))), by=.(type)]
names(temp3) <- c("type", "mean", "se", "lower", "upper")
temp3$election <- "Municipal"
temp <- rbind(temp,temp2,temp3)

quartz(width=5, height=3)
p <- ggplot(data=temp, aes(y=mean, x=type)) + geom_point(aes(color=type)) + theme_bw() + geom_errorbar(aes(ymax=upper, ymin=lower, color=type), width=.2, size=.2) + scale_color_brewer(palette="Dark2", name="Household\ntype") + ylab("Proportion voted\nin any election (2008-2017)") + facet_wrap(~election) + ylim(0,max(temp$upper)+.025) + xlab("Household type") + theme(panel.grid=element_blank())
p

#voting behavior tests
t.test(voteg.prop ~ type, na.rm=TRUE, data=hhdata.matches, paired=TRUE) 
t.test(votep.prop ~ type, na.rm=TRUE, data=hhdata.matches, paired=TRUE)
t.test(votem.prop ~ type, na.rm=TRUE, data=hhdata.matches, paired=TRUE) 



### Figure 6 in paper (voting history by stratum)
temp <- hhdata.matches[,.(mean(voteg.prop, na.rm=TRUE), sd(voteg.prop, na.rm=TRUE)/sqrt(.N), mean(voteg.prop, na.rm=TRUE)-2*(sd(voteg.prop, na.rm=TRUE)/sqrt(.N)), mean(voteg.prop, na.rm=TRUE)+2*(sd(voteg.prop, na.rm=TRUE)/sqrt(.N))), by=.(type, strat)]
names(temp) <- c("type", "strat", "mean", "se", "lower", "upper")
temp$election <- "General"
temp2 <- hhdata.matches[,.(mean(votep.prop, na.rm=TRUE), sd(votep.prop, na.rm=TRUE)/sqrt(.N), mean(votep.prop, na.rm=TRUE)-2*(sd(votep.prop, na.rm=TRUE)/sqrt(.N)), mean(votep.prop, na.rm=TRUE)+2*(sd(votep.prop, na.rm=TRUE)/sqrt(.N))), by=.(type, strat)]
names(temp2) <- c("type", "strat", "mean", "se", "lower", "upper")
temp2$election <- "Primary"
temp3 <- hhdata.matches[,.(mean(votem.prop, na.rm=TRUE), sd(votem.prop, na.rm=TRUE)/sqrt(.N), mean(votem.prop, na.rm=TRUE)-2*(sd(votem.prop, na.rm=TRUE)/sqrt(.N)), mean(votem.prop, na.rm=TRUE)+2*(sd(votem.prop, na.rm=TRUE)/sqrt(.N))), by=.(type, strat)]
names(temp3) <- c("type", "strat","mean", "se", "lower", "upper")
temp3$election <- "Municipal"
temp <- rbind(temp,temp2,temp3)

quartz(width=6, height=3)
p <- ggplot(data=temp, aes(y=mean, x=strat, group=type)) + geom_point(aes(color=type), position=position_dodge(width=0.5)) + theme_bw() + geom_errorbar(aes(ymax=upper, ymin=lower, color=type), width=.2, size=.2, position=position_dodge(width=0.5)) + scale_color_brewer(palette="Dark2", name="Household\ntype") + ylab("Proportion voted\n (2008-2017)") + facet_wrap(~election) + ylim(0,.9) + xlab("Quintile of census tracts by density of solar installations") + theme(panel.grid=element_blank())
p

###


####
# Regressions 
####
library(MuMIn)
library(arm)

## Linear probability model with state fixed effects
hhdata$type.num <- as.numeric(hhdata$type)-1
m.linear <- lm(type.num ~ white.prop + fem.prop + bach.prop + homeown.prop + income_log_hhmean + reg.prop + tsmart_partisan_score_hhmean + voteg.prop + votep.prop + State_name, data=hhdata.matches) #***

stargazer(m.linear, type="text", style="apsr", covariate.labels=c("Prop. White, non-Hispanic", "Prop. female", "Prop. bachelors or higher", "Own home", "Mean HH income (log)", "Prop. registered to vote", "Mean partisan model score", "Prop. voted general election", "Prop. voted primary election"), header=FALSE, dep.var.labels="Household type: solar", star.cutoffs = c(.05, .01, .001),keep=c(1:9))



